#ifndef LBM_UNIT_HPP
#define LBM_UNIT_HPP

#include "LBMUnit.h"

#include <iostream>

namespace ibm
{

template<typename IntType, typename DecType>
LBMUnit<IntType, DecType>::LBMUnit()
{
}
	
template<typename IntType, typename DecType>
LBMUnit<IntType, DecType>::LBMUnit(DecType tau, IntType Nx, IntType Ny, IntType Nz)
:_tau(tau), _Nx(Nx), _Ny(Ny), _Nz(Nz)
{
}
		
template<typename IntType, typename DecType>
LBMUnit<IntType, DecType>::LBMUnit(DecType tau, IntType domain[3], DecType rho_fluid, DecType phy_mu, DecType deltaX)
:_tau(tau), _Nx(domain[0]), _Ny(domain[1]), _Nz(domain[2]), _deltaX(deltaX)
{
	_phy_gravity = -9.81;
}
	
template<typename IntType, typename DecType>
LBMUnit<IntType, DecType>::LBMUnit(DecType tau, std::vector< IntType > domain, DecType rho_fluid, DecType phy_mu, DecType deltaX)
:_tau(tau), _Nx(domain[0]), _Ny(domain[1]), _Nz(domain[2]), _deltaX(deltaX)
{
	_phy_gravity = -9.81;
}
	
template<typename IntType, typename DecType>
LBMUnit<IntType, DecType>::LBMUnit(LBParameterReader<IntType, DecType>* lbparam)
{
	this->_rho_fluid = lbparam->GetFluidDensity();
	this->_phy_mu = lbparam->GetPhysicalDynamicViscosity();
	this->_deltaX = lbparam->GetDeltaX();
	this->_deltaT = lbparam->GetDeltaT();
	this->_tau = lbparam->GetTau();
			
	DecType Lx = lbparam->GetPhysicalLengthX();
	DecType Ly = lbparam->GetPhysicalLengthY();
	DecType Lz = lbparam->GetPhysicalLengthZ();
	
	std::cout << "rho = " << _rho_fluid << std::endl;
	std::cout << "mu = " << _phy_mu << std::endl;
	std::cout << "deltaX = " << _deltaX << std::endl;
	std::cout << "deltaT = " << _deltaT << std::endl;

	this->_L[0] = Lx;
	this->_L[1] = Ly;
	this->_L[2] = Lz;
	
	_phy_gravity = -9.81;




	DecType ReL = (_deltaT / _deltaX) * (1./_deltaX) / (_tau/3. - 1./6.);

	std::cout << "Uf = " << _deltaT / _deltaX << std::endl;
	std::cout << "ReL = " << ReL << std::endl;
	std::cout << "NuL = " << (_tau/3. - 1./6.) << std::endl;


	lbUnits = new olb::LBunits<DecType>(_deltaT / _deltaX, ReL, getResolution(), 1, 1);
}

template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getTau()
{
	//std::cout << lbUnits->getTau() << std::endl;
	//std::cout << _tau << std::endl;
	return this->_tau;
}
	
template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getOmega()
{
	return 1. / this->_tau;
}
	
template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getNu()
{
	//std::cout << lbUnits->getLatticeNu() << std::endl;
	//std::cout << (_tau - 0.5) / 3. << std::endl;
	return (_tau - 0.5) / 3.;
}
	
template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getRe()
{
	std::cout << lbUnits->getRe() << std::endl;
	std::cout << 1. / getNu() << std::endl;
	return 1. / getNu();
}
	
template<typename IntType, typename DecType>
IntType LBMUnit<IntType, DecType>::getNx()
{ 
	return nCell(_L.x()) + 1; 
}
	
template<typename IntType, typename DecType>
IntType LBMUnit<IntType, DecType>::getNy()
{ 
	return nCell(_L.y()) + 1; 
}
	
template<typename IntType, typename DecType>
IntType LBMUnit<IntType, DecType>::getNz()
{ 
	return nCell(_L.z()) + 1; 
}
	
template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getDeltaT()
{
	return this->_deltaT;
}
	
template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getDeltaX()
{
	return this->_deltaX;
}
	
template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getResolution()
{
	return 1. / this->_deltaX;
}
	
template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getFluidDensity()
{
	return this->_rho_fluid;
}
	
template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getPhyNu()
{
	return this->_phy_mu / this->getFluidDensity();
}
	
template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getTimeFactor()
{
	DecType N = getResolution();
	
	DecType nu_ratio = getNu() / getPhyNu();
	
	return N * N / nu_ratio;
	
}
	
template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getLengthFactor()
{
	return getResolution();
}
	
template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getVelocityFactor()
{
	DecType N = getResolution();
	
	DecType nu_ratio = getNu() / getPhyNu();
	
	/*
	std::cout << "N = " << N << std::endl;
	std::cout << "Nu = " << getNu() << std::endl;
	std::cout << "PhyNu = " << getPhyNu() << std::endl;
	*/
	return nu_ratio / N;
}	
	

template<typename IntType, typename DecType>
void LBMUnit<IntType, DecType>::setGravityPhy(DecType grav)
{
	this->_phy_gravity = grav;
}
	
template<typename IntType, typename DecType>
IntType LBMUnit<IntType, DecType>::nCell(DecType L) const 
{ 
	return (int)( L / _deltaX + (DecType)0.5); 
}

template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getGravity()
{
	DecType nu_ratio = getNu() / getPhyNu();
	
	DecType N = getResolution();
	
	DecType gFac = nu_ratio * nu_ratio / (N * N * N);
	
	this->_phy_gravity = -9.81;

	return gFac * _phy_gravity;
}

template<typename IntType, typename DecType>
DecType LBMUnit<IntType, DecType>::getGravityEffected(DecType sphere_density)
{
	DecType nu_ratio = getNu() / getPhyNu();
	
	DecType N = getResolution();
	
	DecType gFac = nu_ratio * nu_ratio / (N * N * N);
	
	DecType density_ratio = getFluidDensity() / sphere_density;
	
	this->_phy_gravity = -9.81;

	return gFac * (1. - density_ratio) * _phy_gravity;
}

} // namespace ibm
#endif // LBM_UNIT_H
